%  Script to read the auroral ground perturbations files send my Dan Weimer
% Feb 24, 2014 Dan Weimer wrote
% Yes, the 360 values are longitudes, starting at 0 and going up to 359
% Yes, there are 90 such rows, going from 0 latitude to 89 latitude.
% I did not do 90 degrees latitude, as all points are the same and N-S and E-W components are zero there.
% And then another header with a time, followed by 90 rows again, for 1440 minutes (records) total.
% A number in the headers that was supposed to be 04 for the month, instead has the minute, so do not use that in this file.


basedirectory = '/Users/manojnair/data/Weimer_Apr_2000/DeltaBresults/unzipped/';

s = dir(basedirectory);

%% East component
ncount = 1;    
for j = 1:10,

    
filename = ['/Users/manojnair/data/Weimer_Apr_2000/DeltaBresults/unzipped/DB_East_200004' sprintf('%02d',j+2) '.txt'];
% East component
fid = fopen(filename,'rt');

for i = 1:1440 , % loop through minute data
    
dummy = fscanf(fid,'%d',5); % get the time header
dummy = fscanf(fid,'%f', 90 * 360);
Y_S(ncount,:,:) = reshape(dummy,[360,90]);

ncount = ncount + 1;
end;

display(filename);
fclose(fid);

end;

save /Users/manojnair/data/Weimer_Apr_2000/DB_Y_20000403_20000412 Y_S sim_time -v7.3; %850 MB

%% North component
ncount = 1;    
for j = 1:10,

    
filename = ['/Users/manojnair/data/Weimer_Apr_2000/DeltaBresults/unzipped/DB_North_200004' sprintf('%02d',j+2) '.txt'];

fid = fopen(filename,'rt');

for i = 1:1440 , % loop through minute data
    
dummy = fscanf(fid,'%d',5); % get the time header
dummy = fscanf(fid,'%f', 90 * 360);
N_S(ncount,:,:) = reshape(dummy,[360,90]);

ncount = ncount + 1;
end;

display(filename);
fclose(fid);

end;

save /Users/manojnair/data/Weimer_Apr_2000/DB_N_20000403_20000412 N_S sim_time -v7.3; %850 MB

%% Vertical component
ncount = 1;    
for j = 1:10,

    
filename = ['/Users/manojnair/data/Weimer_Apr_2000/DeltaBresults/unzipped/DB_Vert_200004' sprintf('%02d',j+2) '.txt'];

fid = fopen(filename,'rt');

for i = 1:1440 , % loop through minute data
    
dummy = fscanf(fid,'%d',5); % get the time header
dummy = fscanf(fid,'%f', 90 * 360);
Z_S(ncount,:,:) = reshape(dummy,[360,90]);

ncount = ncount + 1;
end;

display(filename);
fclose(fid);

end;

save /Users/manojnair/data/Weimer_Apr_2000/DB_Z_20000403_20000412 Z_S sim_time -v7.3; %850 MB

%% plot
for i = 1:10:1440,
imagesc(squeeze(data(i,:,:)));
title(sprintf('%d',i));
pause;

end;
%% Load the observed data and interpolated simulated time series

%load the simulated data
load /Users/manojnair/data/Weimer_Apr_2000/DB_XYZ_20000406

%load obs data
load('/Users/manojnair/projects/GIC/substorm/obs_Apr_2000_all.mat');

% Create a time axis
sim_time = linspace( datenum(2000,4,6,0,0,30), datenum(2000,4,6,23,59,30), 1440);
% name a station

%%
[a1,b1] = sort(theta_gg);

% station = 'NAL';
% [a,b] =intersect(name,station,'rows');


for i = 1:20,
    
 b = b1(i);
% get its location

lon1 = phi_gg(b);
lat1 = 90- theta_gg(b);

% get the magnetic latitude
path(path,'/Users/manojnair/projects/GIC/SHA_obs');
[theta_gm, phi_gm] = gg2gm_new(theta_gg(b),	phi_gg(b));
maglat = - (theta_gm - 90);


% interpolate TS at this observatory (ies)

X_sim = interp3(X_S,ones([1,1440]) * lon1, 1:1440, ones([1,1440]) * lat1 );
Y_sim = interp3(Y_S,ones([1,14400]) * lon1, 1:14400, ones([1,14400]) * lat1 );
Z_sim = interp3(Z_S,ones([1,1440]) * lon1, 1:1440, ones([1,1440]) * lat1 );

subplot(311);
plot(fday, X(b,:) - nanmean(X(b,:)))
hold
plot(sim_time, X_sim,'r');
axis([min(sim_time) max(sim_time) -inf inf]);
datetick('x','keeplimits');
title([name(b,:) sprintf(', MAGLAT = %3.2f',maglat)]);

subplot(312);
plot(fday, Y(b,:) - nanmean(Y(b,:)))
hold
plot(sim_time, Y_sim,'r')
axis([min(sim_time) max(sim_time) -inf inf]);
datetick('x','keeplimits');

subplot(313);
plot(fday, Z(b,:) - nanmean(Z(b,:)))
hold
plot(sim_time, Z_sim,'r')
axis([min(sim_time) max(sim_time) -inf inf]);
datetick('x','keeplimits');

saveas(gcf,['/Users/manojnair/projects/GIC/substorm/Apr2000'  name(b,:) ],'jpeg');

pause;

close;

end;

